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ABSTRACT 

The magnetospheric plasma distribution in the vicinity of a pulsar at various 
inclination angles is investigated using a new relativistic, parallel 3D Particle- 
In-Cell (PIC) code DYMPHNA3D. DYMPHNA3D uses a superposition of the 
electromagnetic fields associated with a rotating magnetised conducting sphere 
in a vacuum (the pulsar fields) and the electromagnetic fields due to the pres- 
ence of the magnetospheric plasma surrounding the pulsar (the plasma fields), 
as the total fields. The plasma is moved self-consistently through the magneto- 
sphere using the PIC methodology Our initial simulation results are presented 
here. These show similar solutions to those obtained from previous numerical 
simulations, which show the fundamental plasma distribution in the vicinity of 
an aligned rotating neutron star to consist of two polar domes and an equa- 
torial torus of trapped non-neutral plasma of opposite sign. The aligned case 
being the case in which the inclination angle, a, between the magnetic dipole 
moment and the rotation axis of the star is zero. Furthermore, our code allows 
for off-axis simulations and we have found that this plasma distribution collapses 
into a 'Quad-Lobe' charge-separated non- neutral magnetospheric plasma distri- 
bution in the case of an orthogonal rotator, i.e., the case in which the magnetic 
dipole moment is oriented at right angles to the rotation axis of the neutron star, 
with the plasma remaining trapped close to the stellar surface by the force-free 
(E • B = 0) surfaces. We find that if initialised with a Goldreich- Julian type dis- 
tribution, the system is seen to collapse rapidly into these stable 'Dome- Torus' 
structures. 

Subject headings: plasmas, methods: numerical, (stars:) pulsars: general 
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Introduction 



The problem as to exactly how and where the observed emission from pulsars is gener- 
ated has defied explanation for nearly four decades. There are many popular models, each 
with its own advantages and disadvantages, each with its dedicated supporters. These mod- 
els differ primarily in terms of the location within the magnetosphere at which the observed 
emission is generated, i.e., the location of the particle accel eration region. The primary high- 
energy emission models at present are the polar c ap model (IHarding fc Muslimov II1998T). the 
outer gap model (ICheng. Ruderman fc Ho Ill986l ). the slot g ap model (IHarding fc Muslimov 
2004 ; IHarding et al. 1120081 ) . and the two-pole caustic model (IDyks fc Rudak 1120031 ). The po- 
lar cap model proposes that the particle acceleration region is located close to the stellar 
surface in the open field region over the pulsar's polar caps. The outer gap model assumes 
that the vacuum gaps within which the particle acceleration occurs form in the outer mag- 
netosphere along null charge surfaces. The slot gap model contains elements of both of these 
models. This model is an extension of the original polar cap model, required to overcome 
some problems that the polar cap model faced in terms of difficulties in producing emission 
beams wide enough to correlate with the wide pulse profiles observed. By allowing the gap to 
extend upwards near the polar cap rim, asymptotically approaching the boundary between 
the open and closed magnetospheric regions, delimited by the last closed field lines, the slot 
gap model is capable of producing these wider pulse profiles. Last, but not least, is the two- 
pole caustic model which assumes that the vacuum gap, in which the particle acceleration 
occurs, is a thin gap that extends from each of the polar caps to the light cylinder and is 
confined to the surface of the last open magnetic field lines. 



Goldreich &: Julian I (119691 ) determined that a magnetised neutron star would be sur- 



rounded 'everywhere' by a charge density given by, 
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where pcj is the Goldreich- Julian charge density and fl is the rotational frequency of the 
star. This was done for the aligned case which is a simplification of the more complex 
inclined rotator problem. They assumed that plasma on open field lines (i.e., the field lines 
which do not close within the light cylinder distance) can stream to infinity, while plasma on 
closed field lines corotates with the neutron star. They also assumed that a wind would be 
generated at the light cylinder, attributable to centrifugal forces due to the rigid corotation 
of this magnetospheric plasma with the star. 

This model advocates a totally filled magnetosphere (TFM), in which the plasma is 
pulled from the pulsar surface to fill the magnetosphere with a charge density given by 
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Eq. (pP). There are some critical problems with this model to include the fact that if the 
plasma has to follow magnetic field lines, then it is impossible for the magnetospheric region 
surrounding an aligned rotator to fill with plasma. Also, at large distance from the stellar 
surface, positive particles would have to travel along magnetic field lines which have their 
feet firmly rooted at the polar regions of the star, where only negative particles exist. The 
assertion that there needs to be acceleration of particles from the polar caps to replace 
particles lost at large distances from the pulsar is unfounded and these issues have been well 
addressed by iMichel I (1l982l . booi l2005h . 



In stark contrast to the scenario in which the pulsar is surrounded 'everywhere' by 
the Goldreich - Julian charge density are t he res ults of explicit computer simulations, first 
performed by iKrause-Polstorff Sz Michel I Jl985al lbl). The results of such simulations demon- 
strate that it is possible to construct stable self-consistent solutions involving charge-separated 
non-neutral plasma distributions in the form of two polar domes of trapped charged parti- 
cles separated by a vacuum gap from an equatorial torus of trapped charged particles of the 
opposite sign, for the case of the aligned rotator, a = 0°. 

The basic premise behind these early simulations was to let the quantised surface charge 
on the neutron star follow the electric forces quasi-statically, using an iterative approach until 
the force-free condition (E • B = 0) is satisfied everywhere. This is an 0(N 2 ) problem, as 
such it is computationally intensive and consequently, the number of particles had to be kept 
reasonably small by choosing a large charge quantisation level. In addition, in order to take 
advantage of the inherent symmetry involved in the system, the charge elements were taken 
to be rings. 



These simulations were later reworked by lSmith et al. I (1200 lh . In this case, the particles 
were taken to be 'super-particles', with each 'super-particle' representing a finite density of 
discrete particles and carrying a charge that is very large in comparison to the electronic 
charge. As before, the code moves the 'super-particles' iteratively to their equilibrium po- 
sitions, i.e., E ■ B = positions based on the electromagnetic fields in the magnetosphere 
which are due to the neutron star and all of the other particles. 

As mentioned, these simulation methodologies are N-body in nature and as a result, 
are restrictive in terms of the number of particles that can be used due to computational 
overheads. The solution in yoked to overcome thi s limita tion is to use a Particle-In-Cell (PIC) 
approach, as was done by ISpitkovsky fc Arons I (120021 ) . To perform their simulations, they 
used a modi fied version of th e publicly available 3D electromagnetic relativistic PIC code 
'TRISTAN' (IBuneman Ill993l ). Their approach used a superposition of the vacuum rotator 
fields, i.e., th e electromagnet ic fields associated with a rotating conducting magnetised sphere 
in a vacuum (IDeutsch Ill955l ) and the plasma fields, i.e., the electromagnetic fields due to the 
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presence of the magnetospheric plasma, as the total electromagnetic fields. Deutsch wrote 
down by inspection, the vacuum rotator field equations (providing little insight into their 
derivation) which require a significant amount of processing to convert them i nto functions 
of r, 9 & 4>. An explicit derivation of these field equations was performed by iMichel fc Li 
f)l999h . Eq. ©. 
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where B is the magnetic field strength at the equator, r is the radial distance from the 
star, 9 is the stellar latitude, (ft is the stellar longitude (</> s = — Qt), and a is the inclination 
angle between the magnetic dipole moment and the rotation axis of the neutron star. 

Our simulations also use 'super-particles', each of which represents a finite density of 
discrete particles and carries a charge that is very large in comparison to the electronic 
charge. The aligned rotator scenario mentioned above is the case of maximal symmetry. It 
is a simplification of the significantly more complicated inclined rotator scenario, in which the 
magnetic dipole moment is oriented at an angle, a, to the rotat ion axis, i.e., an orthogonal 



rotator being the case of minimal symmetry. It has been argued ( IMichel 1120051 ) that focusing 
efforts on the aligned case is a futile endeavour as this case does not produce an active pulsar. 
Also, the physics of this scenario is so drastically different to that of the inclined/orthogonal 
case, that the aligned model is incapable of generating any real insight into the operation of 
pulsars. Instead it is contested that the inclination between the magnetic a nd rotational axes 



of the neutron star is an essential attribute of any (active) pulsar model (IMichel fc Smith 
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200ll ; Michel 1120051 ) 



These simulations have independently verified the original KPM simulation results and 
showed the magnetospheric charge distribution to consist of these two polar domes of trapped 
negative charge (electrons), with an equatorial torus of trapped positive charge (protons, 
positrons or ions), separated by a vacuum gap. It is worth noting that the fundamentals of 
this 'Dome- Torus' ma gnetospheric plasma distribution hav e also been i nvesti g ated and ver 



ified in dependently bvlRvlov I Jl976l. Il989h ; Ishibata I Jl989[) 



f|1993h 
J2004T ) 



Thielheim fc Wolfsteller 



(11994h . and IPetri et al. 



Neukirch I (119931 ); IZachariades 



I2OO2J) as indicated by |Mchel 



2. Simulation code 

The code currently under development is a 3D fully electromagnetic, parallel, scalable, 
and relativistic PIC code, which has been designed in an extremely modular fashion. The 
codes modular design will facilitate ease of application of the code to a wide variety of kinetic 
plasma phenomena with little modification. 

Due to its inherent simplicity, PIC methodology is widely used for the kinetic simulation 
of plasmas. Its implementation involves the imposition of a spatial mesh (a spatial discreti- 
sation) on the reg ion containing the plasma. A staggered mesh system known as the Yee 



lattice (lYee Ill966l ) is commonly utilised for this purpose. Electromagnetic field quantities 
are defined on the spatial grid whereas particles can occupy arbitrary locations in contin- 
uum space. The attributes of the plasma, i.e., the charge density and current density (field 
sources) are deposited onto the spatial mesh and thus, the plasma particles need only interact 
with the mesh, as opposed to each other as in the case of a Coulomb problem. This reduces 
a computationally infeasible 0(N 2 ) problem to a computationally viable 0(N log N) scale. 
Once the plasma attributes have been deposited onto the spatial mesh, Maxwell's equations 
are self-consistently solved on the mesh to determine the electromagnetic fields due to the 
presence of the plasma. These plasma fields are then superimposed with the vacuum rotator 
(pulsar) fields, i.e., the Deutsch fields, to determine the total electromagnetic fields. The 
total fields, located on the mesh, are then interpolated back to the location of the particles, 
the Lorentz force calculated and the particles are moved by integration of Newton's laws 
of motion. This procedure is iterated (a temporal discretisation) self-consistently moving 
the plasma through the computational domain in accordance with the total fields calculated 
above. 
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2.1. Data output 

One of the major difficulties with modern large-scale simulation codes is the sheer volume 
of the output data produced in many cases. DYMPHNA3D is no exception and is capable 
of producing a very large amount of output data, depending on the scale of the simulation 
and the configuration of the desired output. With this in mind, throughout the development 
and construction of the code, we focused heavily on the data output for post-processing 
and visualisation of the results in order to take maximal advantage of the 3D nature of the 
simulation and to be able to efficiently process the output. 

The user has the option of outputting the required simulation data at the end of each 
iteration (or specified increment of iterations), in order to study the temporal evolution or 
just the final state at the end of the simulation run if that is all that is required. The choice 
of data to output includes the Cartesian components of the particle locations, velocities, the 
Lorentz force on the particles, the electric field, the magnetic field, and the current density. 
The particle energies, their gamma factors, E • B on the mesh, and the charge density on 
the grid as well as specific phase space data for the particles can also be output if required. 
All data can be output in either standard data files for arbitrary post-processing or in VTK 



(Visualisation Toolkit, http://www.vtk.org) file format. 



3. Numerical model 

As an initial test of the code, the investigation of the magnetospheric plasma distribution 
in the vicinity of a pulsar was chosen. The PIC code has been imposed on the region 
surrounding the pulsar, having the star centred in the computational domain. For simplicity, 
we have assumed that the positive and negative charge carrying species are of equal mass, 
alleviating any mass disparity complications in the computational model. We have initialised 
a neutral plasma (surface charge) on the pulsar surface, i.e., positive and negative charge 
carrying 'super-particle' pairs are initialised at the same location on the stellar surface. This 
is done in order to satisfy initial conditions requiring charge neutrality, such that Gauss' law, 
Eq. 0, is satisfied initially and as long as the charge continuity equation, Eq. (TJJ, is not 
violated, it will be satisfied at subsequent times during the simulation. 

V-E = p (3) 



(4) 
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Taking Maxwell's V-B = equation as a second initial condition, leaves Maxwell's two curl 
equations, from which a set of update equations for the electromagnetic fields due to the 
presence of the magnetospheric plasma can be derived, Eqs. IA1HA61 



3.1. Simulations 

Simulations were run on the 32 processor BULL Novascale shared-memory machine 
at the ICHEC 1 facility. As discussed above, the total electromagnetic fields utilised in the 
simulations are a superposition of the plasma fields and the pulsar fields. The electromagnetic 
field structure in the vicinity of the pulsar has been visualised by means of 3D vector plots 
using Paraview ( http: / / www.paraview.org]) . Figures [if & [2f illustrate the structure of the 



fields for various angles of inclination, a, ranging from the aligned rotator scenario, a = 0° 
(left), to the case of the orthogonal rotator, a = 90° (right). The images in these two figures 
are taken looking in along the negative y-axis towards the origin. The dipolar nature of the 
magnetic field is overtly evident in the leftmost panel of Fig. ([T]) and is seen to retain this 
dipolar form throughout the inclination from aligned rotator, a = 0°, to orthogonal rotator, 
a = 90°. 

The case however, is very different for the electric field during inclination. In the leftmost 
panel of Fig. (J2J), i.e., for the aligned rotator scenario, the structure of the E field can be 
observed to be primarily quadrupolar in nature, but its structure is seen to evolve significantly 
with inclination. From a visual inspection of the set of equations in Eq. (T5]), it can clearly 
be seen that they are simply a superposition of the aligned (cos a terms, symmetrical field 
structure, with no time-dependence) and orthogonal (sin a terms, anti-symmetrical field 
structure, with time-dependence) cases. The simplicity of the aligned rotator scenario with 
respect to that of the orthogonal rotator is acutely evident from such an analysis of the field 
equations, illustrating the reasoning behind the historically predominant investigation of the 
former. 

In terms of our investigation of the plasma distribution in the vicinity of a pulsar, 
the initialised neutral plasma (surface charge) on the stellar surface is released into the 
magnetosphere and it self-consistently follows the ambient fields driven by the iterative 



1 Irish Centre for High End Computing, http: / / www.ichecle] 



2 Note that it is the length of the field vectors and their respective directions that gives the illusion of 
anti-symmetry in the images, the fields are symmetrical. We plot only a portion of the total number of field 
vectors, for the purpose of clarity, which also aids in the illusion of anti-symmetry. The purple axis denotes 
the stellar rotation axis while the white axis marks the magnetic dipole moment. 
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PIC algorithm in accordance with the total determined fields, as discussed above. Once 
released, the particles are numerically free to travel anywhere within the computational 
domain (governed by the ambient electromagnetic fields), which includes reentering the star 
if they want to. With typical simulation runs containing of the order of 2M to 20M 'super- 
particles', there is always a supply of charge available on the stellar surface for emission into 
the magnetosphere, wherever it is needed by the ambient fields. As mentioned above, the 
code can output the state of the simulation at each iteration (or increment of iterations), 
which allows us to study the temporal and spatial evolution of the magnetospheric plasma 
distribution. 

The simulation is terminated when a stable configuration for the non-neutral plasma is 
obtained. The condition for this is that the plasma has distributed itself, such that it elim- 
inates the parallel component of the electric field along the magnetic field lines responsible 
for pulling the charges from the stellar surface and accelerating them out along the mag- 
netic field lines. Thus, some form of trigger is required to terminate the simulation when 
this condition has been reached numerically. There are two relatively straight forward ways 
of testing for the presence of this stable configuration, firstly one can monitor the value of 
E • B inside the plasma itself, which should be zero everywhere within the plasma for it to 
have equilibriated. Secondly, the parallel component of the electric field along the magnetic 
field lines serves to pull charges from the stellar surface as mentioned above. Thus, if this 
parallel component of the electric field has been eliminated by the plasma redistributing 
itself to do so, then the existence of this condition will be marked by a cessation of particle 
emission from the st e llar s urface. This was essentially, the criterion that was utilised by 



Spitkovsky fc Arons I (120021 ). to identify when this configuration had been achieved in their 



simulations and is also the trigger that we have employed in our model. 

As the simulation progresses, the code continually monitors the motion of particles 
released from the stellar surface for a cessation of emission into the magnetosphere. Once the 
code has determined that a cessation of emission of plasma has been achieved, i.e., particles 
are no longer being pulled from the surface of the star as a result of the parallel component 
of the electric field above the surface being eliminated, this triggers the termination of the 
simulation. We find that the stable equilibrium in the case of an aligned rotator, a = 0°, is in 
the form of a charge-separated distribution consisting of two polar domes and an equatorial 
torus of trapped non-neutral plasma of opposite signs. 
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3.2. Force-Free surfaces - Trapping regions 

The fundamental reason that we see these stable charge-separated non-neutral plasma 
structures is clearly evident when we look at E ■ B in the vicinity of the neutron star, taking 
a slice through the pulsar origin in the x-z plane for clarity. Plots of E • B in the vicinity of 
the stellar surface for various inclination angles, a, have been included in Fig. Q. It should 
be noted that the images of E • B in Fig. (j3J) depict the vacuum rotator fields only, i.e., in the 
absence of plasma surrounding the star. In these images, red indicates that E • B > 0, blue 
indicates that E • B < and green indicates E • B = 0. The geometry o f these force-free 



surfaces is in good agreement with those illustrated in iMichel fc Li I (119991 ). It can be seen 
in the aligned case, that the polar domes are a reflection of these force-free trapping surfaces 
over the caps of the pulsar with negative charge being extracted from the poles, driven by 
the quadrupolar electric field. In the case of the equatorial torus, there is a force-free surface 
along the equatorial plane, on which the particles would tend to gather if it were not for their 
inter-repulsion which pushes them out to form the observed torus constrained outwardly by 
the tightly curved magnetic field lines in this region. 

In the case of the orthogonal rotator, depicted in the rightmost panel of Fig. ([3]), the 
'Quad-Lobe' structure is clearly evident in the locus of these force-free trapping regions in 
the vicinity of the stellar surface. 



3.3. Plasma distribution - Simulation results & Discussion 

An illustration of the final particle locations predicted from our simulations has been 
included in Fig. (jlj), for the same angles of inclination included for the previous images, 
a = 0°, 45° and 90°. In the case of the aligned rotator, the 'Dome- Torus' charge-separated 
plasma distribution has been found to be the stable equilibrium condition, in agreement with 
the results of the previous numerical simulations performed. 

As discussed above, this 'Dome- Torus' distribution is not seen to persist with rotation 
of the magnetic dipole moment with respect to the stellar rotation axis. In the case of 
inclination to the orthogonal rotator scenario, i.e., a = 90°, the 'Dome- Torus' structure is 
observed to collapse into a 'Quad-Lobe' type distribution as illustrated in the rightmost panel 
of Fig. (j4j). This 'Quad-Lo be' structure bears a stark resemblance to the plasma distribution 



predicted by lParish I (119741 ) (charge density given by Eq. (jSJ)), with the exception that instead 
of being allowed to propagate out to a significant percentage of the light cylinder distance 
as was predicted in that early theoretical work, the plasma remains bound much closer to 
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the stellar surface by these force-free trapping surfaces. 



P 
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(5) 



In the intermediate case, in which the inclination angle a = 45°, the negative charge 
carrying species is seen to occupy the lopsided dome evident in the locus of the polar force- 
free surfaces at this inclination in Fig. Q. The positive charge carrying species is seen 
to propagate to a greater radial distance along the vertically displaced E • B = surfaces 
that lie along the equatorial plane in the aligned case and are found to be stable in this 
configuration. 

We have also found that if initialised with the Goldreich- Julian TFM state, the system 
is unstable and collapses into these stable trapped non-neutral plasma structures as was first 
demonstrated numerically by ISmith et al. I (120011 ). 



The fundamental problem with this scenario, in which the plasma is trapped in these 
non-neutral structures close to the stellar surface, is that observational evidence clearly 
illustrates that in practise, the plasma is not confined to these trapping zones. Images of 
the Crab nebula unambiguously depict a large-scale equatorial toroidal ou tflow and polar 
jets reminiscent to the topology of the 'Dome- Torus' structure shown here. iMichel I (120051 ) 
theorises that wave pickup of the charged particles by the large amplitude electromagnetic 
waves generated by the orthogonal component of the magnetic dipole field could be the 
resolution to this particular issue. 



4. Conclusions 

The initial test application of this newly developed code has yielded some interesting 
results. We have found from our investigations into the plasma distribution in the vicinity 
of a rapidly rotating neutron star, that stable solutions involving charge-separated non- 
neutral plasma structures are not only viable, but are the natural equilibrium states for 
these systems. Our results agree well with the results of prior numerical simulations based 
on the aligned rotator scenario and we have extended these to look at the 'off-axis' cases as 
well. 

The general topology of the Crab nebula, as observed in images from the Chandra X-ray 
observatory, is in some fashion reminiscent of the 'Dome- Torus' type plasma distribution, 
evident in the observed equatorial toroidal outflow and polar jets. If the observed equatorial 
toroidal outflow and polar jets are in some way a reflection of this kind of plasma structure 
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in the vicinity of the pulsar, then given that this distribution is seen to collapse into a 
'Quad-Lobe' type structure in the case of the orthogonal rotator, this may suggest that a 
nebula involving an orthogonal rotator or near orthogonal rotator, would not posses this 
observable macroscopic structure. This would also indicate that the inclination of the Crab 
pulsar should be appreciably lower than than that of an orthogonal rotator, in order that 
it still possess this macroscopic quasi-'Dome- Torus' structure. Chandra observations of the 
young Crab-like pulsar J0205+6449 in the SNR 3C 58 also clearly show a macroscopic torus 
structure (seen edge-o n) and some j et-lik e stru cture. The inclination angle is estimated to 
be of the order of 70°, Slane et al. (|2002 . 2004 ). slightly in excess of typical Crab estimates 
of around 65°, Izhang fc Cheng I J2002h . 



Given the success of this fledgling code in its initial application and its inherent mod- 
ularity, these pave the way for application to a plethora of interesting astrophysical kinetic 
plasma phenomena. Future plans for the code involve continuing the initial investigations 
presented in this pape r in greater detail to include probing the possibility and effects of the 
diocotron instability, (jSpitkovsky fc Arons 1 120021 ; iPetri et al. 1 120021 . 120031 ). Also planned as 
part of this work is the investigation of possible mechanisms required to perturb the observed 
stable distribution in order to allow the plasma to propagate beyond these trapping regions 
to form the observed equatorial outflow and polar jets. The nature of the developed code 
will allow us to numerically investigate the pulsar magnetospheric plasma distribution in 
unprecedented detail both spatially and temporally. Investigations of possible pulsar GRP 
mechanisms, necessarily over very short spatial scales are planned as well as simulations of 
a mechanism capable of generating broadband persistent electron cyclotron maser emission 
from Brown Dwarfs. The developed code is at its most basic, a 'first-principles' approach to 
modelling plasmas and has been implemented in an extremely modular fashion facilitating 
the possibility of inclusion a wide variety kinetic phenomena into the computational model. 
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A. Plasma fields update equations 
The update equation for E x : 
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The update equation for E y : 
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The update equation for E z : 
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The update equation for B x : 
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The update equation for B y : 
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The update equation for B z : 
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Fig. 1. — 3D vector plots of the magnetic field structure in the vicinity of a pulsar for 
various angles of inclination, a, from Paraview, ranging from the aligned case (leftmost 
panel) through 45° inclination (centre panel) to the orthogonal case (rightmost panel). The 
colour scheme tends from red (maximum magnitude vector facing upwards) through orange, 
yellow, centred on green, through turquoise and light blue to dark blue (maximum magnitude 
vector facing downwards). It is evident from these images that the magnetic field is seen to 
retain its dipolar form throughout the inclination from aligned to orthogonal rotator. Note 
that it is the length of the field vectors and their respective directions that gives the illusion 
of anti-symmetry in the images, the fields are symmetrical. We plot only a portion of the 
total number of field vectors, for the purpose of clarity, which also aids in the illusion of 
anti-symmetry. The purple axis denotes the stellar rotation axis while the white axis marks 
the magnetic dipole moment. Image taken looking in along the negative y-axis towards the 
origin. 



Fig. 2. — 3D vector plots of the electric field structure in the vicinity of a pulsar for vari- 
ous angles of inclination, a, from Paraview, ranging from the aligned case (leftmost panel) 
through 45° inclination (centre panel) to the orthogonal case (rightmost panel). The colour 
scheme tends from red (maximum magnitude vector facing upwards) through orange, yellow, 
centred on green, through turquoise and light blue to dark blue (maximum magnitude vector 
facing downwards). It is clear from the images that there is a significant evolution of the 
electric field structure from its initial quadrupolar nature during the inclination. We plot 
only a portion of the total number of field vectors, for the purpose of clarity. The purple axis 
denotes the stellar rotation axis while the white axis marks the magnetic dipole moment. 
Image taken looking in along the negative y-axis towards the origin. 
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Fig. 3. — Plots of E • B in the vicinity of a pulsar for various inclination angles, a, from 
Paraview, ranging from the aligned case (leftmost panel) through 45° inclination (centre 
panel) to the orthogonal case (rightmost panel). Red indicates that E • B > 0, blue 
indicates that E ■ B < and green indicates E ■ B = 0. Images are taken in a slice through 
the pulsar origin in the x-z plane. The shape of the trapping regions is evident in the locus 
of the E ■ B = surfaces in these plots. The purple axis denotes the stellar rotation axis 
while the white axis marks the magnetic dipole moment. It should be noted that this image 
depicts the vacuum rotator fields only, i.e., in the absence of plasma surrounding the star. 
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Fig. 4. — Plots of final particle locations in the vicinity of a pulsar for various inclination 
angles, a, from Paraview, ranging from the aligned case (leftmost panel) through 45° incli- 
nation (centre panel) to the orthogonal case (rightmost panel). Red indicates the negative 
charge carrying species and yellow indicates the positive charge carrying species. The 'Dome- 
Torus' structure can be seen in the aligned case in the leftmost panel and can be seen to 
collapse into this 'quad-lode'-like structure in the orthogonal case in the rightmost panel. 
The purple axis denotes the stellar rotation axis while the white axis marks the magnetic 
dipole moment. Image taken looking in along the negative y-axis towards the origin. 



